# show estimation results

using FileIO, JLD2

include("cenf.jl")

cenf.readData(root)

f=load("$(root)/output/ms_f_output_withfixtheta_refined.jld2")
ms_fvals = f["ms_fvals"]
ms_θ = f["ms_θ"]
ms_β = f["ms_β"]
ms_η = f["ms_η"]
ms_γ = f["ms_γ"]
ms_logT = f["ms_logT"]
ms_logS = f["ms_logS"]

idx_best = findfirst(isequal(maximum(ms_fvals)), vec(ms_fvals));

bestθ = ms_θ[idx_best]
bestβ = vec(ms_β[idx_best,:])
bestη = vec(ms_η[idx_best,:])
bestγ = reshape(ms_γ[idx_best,:],35,35)
bestlogT = reshape(ms_logT[idx_best,:],109,35)
bestlogS = reshape(ms_logS[idx_best,:],109,35)

oP=cenf.solveWithFixedTheta(2.0,bestβ,bestη,bestγ,bestlogT,bestlogS,"fitted_FF")
println("------------------------------------")
println("Fixed θ, Z_F")
println("beta[1]: ", oP["β"][1])
println("         (", sqrt(oP["Var_β1"]),")")
println("beta[2]: ", oP["β"][2])
println("         (", sqrt(oP["Var_β2"]),")")
println("eta[1]: ", oP["η"][1])
println("        (", sqrt(oP["Var_η1"]),")")
println("eta[2]: ", oP["η"][2])
println("        (", sqrt(oP["Var_η2"]),")")
println("pseudo-loglikelihood: ", -oP["fval"])
println("pseudo-r2           : ", oP["pseudor2"])
println("------------------------------------")
oPFF = oP


include("cenf.jl")
cenf.readData(root)

f=load("$(root)/output/ms_f_output_withvartheta_refined.jld2")
ms_fvals = f["ms_fvals"]
ms_θ = f["ms_θ"]
ms_β = f["ms_β"]
ms_η = f["ms_η"]
ms_γ = f["ms_γ"]
ms_logT = f["ms_logT"]
ms_logS = f["ms_logS"]

idx_best = findfirst(isequal(maximum(ms_fvals)), vec(ms_fvals));

bestθ = [ms_θ[idx_best]]
bestβ = vec(ms_β[idx_best,:])
bestη = vec(ms_η[idx_best,:])
bestγ = reshape(ms_γ[idx_best,:],35,35)
bestlogT = reshape(ms_logT[idx_best,:],109,35)
bestlogS = reshape(ms_logS[idx_best,:],109,35)

oP=cenf.solveWithθ_f(bestθ,bestη,bestγ,bestlogT,bestlogS,"fitted_VF")
println("------------------------------------")
println("Variable θ, Z_F")
println("theta  : ", oP["θ"])
println("         (", sqrt(oP["Varθ"]),")")
println("eta[1]: ", oP["η"][1])
println("        (", sqrt(oP["Var_η1"]),")")
println("eta[2]: ", oP["η"][2])
println("        (", sqrt(oP["Var_η2"]),")")
println("pseudo-loglikelihood: ", -oP["fval"])
println("pseudo-r2           : ", oP["pseudor2"])
println("------------------------------------")
oPVF = oP

using Printf
# produce core of a latex table that shows the results
open("$(root)/output/structestout.tex", "w") do fh
  write(fh, @sprintf("\$\\theta\$ & 4 & %0.2f \\\\ \n",oPVF["θ"]))
  write(fh, @sprintf("& & (%0.2f) \\\\ \n",sqrt(oPVF["Varθ"])))
  write(fh, @sprintf("\\midrule \$z\$ measure used & \$z^{(2)}\$ & \$z^{(2)}\$  \\\\ \n"))
  write(fh, @sprintf("\\midrule Pseudo-loglikelihood & %0.2f & %0.2f \\\\ \n", oPFF["fval"]  , oPVF["fval"]))
  write(fh, @sprintf("Pseudo-\$R^2\$ & %0.2f & %0.2f \\\\ \n",oPFF["pseudor2"] ,oPVF["pseudor2"]))
end

# produce core of table
using Statistics
open("$(root)/output/structestout_detail.tex", "w") do fh
  write(fh, @sprintf("""
\$ \\beta_1\$ & %0.2f   & \\\\
		      & (%0.2f)  & \\\\
\$ \\beta_2\$ & %0.2f    & \\\\
		      & (%0.2f)  & \\\\
\$ \\eta_1\$  & %0.2f  & %0.2f\\\\
		      & (%0.2f) & (%0.2f) \\\\
\$ \\eta_2\$  & %0.4f & %0.4f \\\\
          & (%0.4f) & (%0.4f)  \\\\
""",  oPFF["β"][1], 
      sqrt(oPFF["Var_β1"]), 
      oPFF["β"][2], 
      sqrt(oPFF["Var_β2"]) , 
      oPFF["η"][1]*std(vec(cenf.z_f[1,:,:])), oPVF["η"][1]*std(vec(cenf.z_f[1,:,:])) , 
      sqrt(oPFF["Var_η1"])*std(vec(cenf.z_f[1,:,:])), sqrt(oPVF["Var_η1"])*std(vec(cenf.z_f[1,:,:])), 
      oPFF["η"][2]*std(vec(cenf.z_f[1,:,:]))^2, oPVF["η"][2]*std(vec(cenf.z_f[1,:,:]))^2,
      sqrt(oPFF["Var_η2"])*std(vec(cenf.z_f[1,:,:]))^2, sqrt(oPVF["Var_η2"])*std(vec(cenf.z_f[1,:,:]))^2))
end
